#define SOL_ALL_SAFETIES_ON 1
#include <sol/sol.hpp>
#include "matrix.hpp"
#include "vector.hpp"
#include "tensor.hpp"

#include <iostream>
#include <list>
#include <cmath>
#include <complex>
#include "complex.h"
#include <omp.h>
#include <log/log.h>
#include <appendix/phy_constants.h>
#include <boost/math/special_functions/sin_pi.hpp>
#include <boost/math/special_functions/sinc.hpp>
namespace sol {
   template <>
   struct is_automagical<std::vector<std::complex<double> > > : std::false_type {};
}

typedef Eigen::Tensor<double,3> Tensor3;
typedef Eigen::Tensor<double,4> Tensor4;
typedef Eigen::Tensor< std::complex<double> ,4> CTensor4;
typedef Eigen::Tensor< std::complex<double> ,3> CTensor3;
typedef Eigen::Matrix<std::string, Eigen::Dynamic, Eigen::Dynamic> EigenSMatrix;

namespace sol {
        //template <>
        //struct is_automagical<EigenMatrix> : std::false_type {};
        template <>
        struct is_automagical<Tensor4> : std::false_type {};
        template <>
        struct is_automagical<Tensor3> : std::false_type {};
        template <>
        struct is_automagical<CTensor4> : std::false_type {};
        template <>
        struct is_automagical<CTensor3> : std::false_type {};
}

namespace LuaNumeric {
	typedef std::complex<double> complex;
	typedef std::vector<double> vector;
	typedef std::vector<int> ivector;
	typedef std::vector<complex> cvector;
	typedef std::vector<vector> vvector;
	typedef std::vector<ivector> vivector;

	void open_luaGlobal(sol::table& module) {
		module.set_function("num_threads", sol::overload(
				[](int num){
					logs::info_log(boost::format("set threads: %1%") % num);
					omp_set_num_threads(num);
					return omp_get_num_threads();
				},
				[](){
					return omp_get_num_threads();
				}
			)
		);

		module.set_function("grids", sol::overload(
								[](const vector& x, const vector& y)
								{
									auto xx = std::make_shared<vector>( x.size()*y.size() );
									auto yy = std::make_shared<vector>( x.size()*y.size() );
									int xnum = x.size();
									int ynum = y.size();

									for( int ii=0; ii < xnum; ii ++)
									{
										for( int jj=0; jj < ynum; jj ++)
										{
											(*xx)[ii*ynum+jj] = x[ii];
											(*yy)[ii*ynum+jj] = y[jj];
										}
									}
									return std::make_tuple(xx,yy);
								},
								[](const cvector& x, const cvector& y)
								{
									auto xx = std::make_shared<cvector>( x.size()*y.size() );
									auto yy = std::make_shared<cvector>( x.size()*y.size() );
									int xnum = x.size();
									int ynum = y.size();

									for( int ii=0; ii < xnum; ii ++)
									{
										for( int jj=0; jj < ynum; jj ++)
										{
											(*xx)[ii*ynum+jj] = x[ii];
											(*yy)[ii*ynum+jj] = y[jj];
										}
									}
									return std::make_tuple(xx,yy);
								},
								[](const vvector& vec,vvector& g, vivector& idx)
								{
									g.clear();
									idx.clear();
									Vector::grids_index(vec, idx);
									g.reserve(idx.size());

									for(auto&& id:idx){
										vector d(id.size());
										for(int jj=0;jj<id.size(); jj ++ )
										{
											d[jj] = vec[jj][id[jj]];
										}
										g.push_back(d);
									}
								},
								[](const vvector& vec)
								{
									vvector g;
									vivector idx;
									Vector::grids_index(vec, idx);
									g.reserve(idx.size());
									for(auto&& id:idx){
										vector d(id.size());
										for(int jj=0;jj<id.size(); jj ++ )
										{
											d[jj] = vec[jj][id[jj]];
										}
										g.push_back(d);
									}
									return std::make_tuple(g,idx);
								}
			)
		);
		module.set_function("serialize", sol::overload(
					[]( const std::vector<double>& vec ){
						return Vector::serialize(vec);
					},[]( const std::vector<std::complex<double> >& vec ){
						return Vector::serialize(vec);
					},[]( const std::complex<double>& v ){
						return serialize(v);
					},[]( const Tensor3& v ){
						return Tensor::serialize(v);
					},[]( const Tensor4& v ){
						return Tensor::serialize(v);
					},[]( const CTensor3& v ){
						return Tensor::serialize(v);
					},[]( const CTensor4& v ){
						return Tensor::serialize(v);
					},[]( const EigenSMatrix& m ){
						return Eigen::serialize(m);
					}
				)
		);

		module.set_function("grids_idx", sol::overload(
						[](const vvector& vec)
						{
							vivector idx;
							Vector::grids_index(vec, idx);
							return idx;
						},
						[](const vvector& vec, vivector& res)
						{
							Vector::grids_index(vec, res);
						}
				)
		);


		module.set_function("kronecker_delta", [](int ii,int jj){ return double(ii==jj); } );

		module.set_function("round", [](double x){ return int( std::round(x) ); } );
		module.set_function("abs", sol::overload(
						sol::resolve<double(double)>(std::abs),
						sol::resolve<double(const complex&)>(std::abs<double>),
						sol::resolve<vector(const cvector&)>(Vector::abs<complex>),
						sol::resolve<vector(const vector&)>(Vector::abs<double>) )
		);

		module.set_function("sqrt", sol::overload(
						sol::resolve<double(double)>(std::sqrt),
						sol::resolve<complex(const complex&)>(std::sqrt<double>),
						sol::resolve<cvector(const cvector&)>(Vector::sqrt<complex>),
						sol::resolve<vector(const vector&)>(Vector::sqrt<double>) )
		);

		module.set_function("dot", sol::overload(
						sol::resolve<double(const vector&,const vector&)>(Vector::dot<double>),
						sol::resolve<complex(const cvector&,const cvector&)>(Vector::dot<complex>)
						)
		);
		module.set_function("sum", sol::overload(
						sol::resolve<double(const vector&)>(Vector::sum<double>),
						sol::resolve<complex(const cvector&)>(Vector::sum<complex>)
						)
		);
		module.set_function("real", sol::overload(
						sol::resolve<double(const complex&)>(std::real<double>),
						sol::resolve<vector(const cvector&)>(Vector::real<complex>),
						sol::resolve<vector(const vector&)>(Vector::real<double>)
						)
		);
		module.set_function("conj", sol::overload(
						sol::resolve<complex(const complex&)>(std::conj<double>),
						sol::resolve<cvector(const cvector&)>(Vector::conj<complex>)
						)
		);
		module.set_function("imag", sol::overload(
						sol::resolve<double(const complex&)>(std::imag<double>),
						sol::resolve<vector(const cvector&)>(Vector::imag<complex>),
						sol::resolve<vector(const vector&)>(Vector::imag<double>)
						)
		);
		module.set_function("arg", sol::overload(
						sol::resolve<double(const complex&)>(std::arg<double>),
						sol::resolve<vector(const cvector&)>(Vector::arg<complex>),
						sol::resolve<vector(const vector&)>(Vector::arg<double>),
						sol::resolve<vector(const vector&,const vector&)>(Vector::arg)
						)
		);
		module.set_function("phase", sol::overload(
						sol::resolve<vector(const cvector&)>(Vector::phase<complex>) )
		);

		module.set_function("log", sol::overload(
						sol::resolve<double(double)>(std::log),
						sol::resolve<complex(const complex&)>(std::log<double>),
						sol::resolve<cvector(const cvector&)>(Vector::log<complex>),
						sol::resolve<vector(const vector&)>(Vector::log<double>) )
		);

		module.set_function("log10", sol::overload(
						sol::resolve<double(double)>(std::log10),
						sol::resolve<complex(const complex&)>(std::log10<double>),
						sol::resolve<cvector(const cvector&)>(Vector::log10<complex>),
						sol::resolve<vector(const vector&)>(Vector::log10<double>) )
		);
		module.set_function("sin", sol::overload(
						sol::resolve<double(double)>(std::sin),
						sol::resolve<complex(const complex&)>(std::sin<double>),
						sol::resolve<cvector(const cvector&)>(Vector::sin<complex>),
						sol::resolve<vector(const vector&)>(Vector::sin<double>) )
		);
		module.set_function("sinc_pi", sol::overload(
						sol::resolve<vector(const vector&)>(Vector::sinc_pi),
						//sol::resolve<cvector(const cvector&)>(Vector::sinc_pi),
						sol::resolve<double(double)>(boost::math::sinc_pi<double>)
		));
		module.set_function("cos", sol::overload(
					sol::resolve<double(double)>(std::cos),
						sol::resolve<complex(const complex&)>(std::cos<double>),
						sol::resolve<cvector(const cvector&)>(Vector::cos<complex>),
						sol::resolve<vector(const vector&)>(Vector::cos<double>) )
		);
		module.set_function("asin", sol::overload(
						sol::resolve<double(double)>(std::asin),
						sol::resolve<complex(const complex&)>(std::asin<double>),
						sol::resolve<cvector(const cvector&)>(Vector::asin<complex>),
						sol::resolve<vector(const vector&)>(Vector::asin<double>) )
		);
		module.set_function("acos", sol::overload(
						sol::resolve<double(double)>(std::acos),
						sol::resolve<complex(const complex&)>(std::acos<double>),
						sol::resolve<cvector(const cvector&)>(Vector::acos<complex>),
						sol::resolve<vector(const vector&)>(Vector::acos<double>) )
		);

		module.set_function("tan", sol::overload(
						sol::resolve<double(double)>(std::tan),
						sol::resolve<complex(const complex&)>(std::tan<double>),
						sol::resolve<cvector(const cvector&)>(Vector::tan<complex>),
						sol::resolve<vector(const vector&)>(Vector::tan<double>) )
		);

		module.set_function("atan", sol::overload(
						sol::resolve<double(double)>(std::atan),
						sol::resolve<complex(const complex&)>(std::atan<double>),
						sol::resolve<cvector(const cvector&)>(Vector::atan<complex>),
						sol::resolve<vector(const vector&)>(Vector::atan<double>) )
		);

		module.set_function("sinh", sol::overload(
						sol::resolve<double(double)>(std::sinh),
						sol::resolve<complex(const complex&)>(std::sinh<double>),
						sol::resolve<cvector(const cvector&)>(Vector::sinh<complex>),
						sol::resolve<vector(const vector&)>(Vector::sinh<double>) )
		);
		module.set_function("cosh", sol::overload(
						sol::resolve<double(double)>(std::cosh),
						sol::resolve<complex(const complex&)>(std::cosh<double>),
						sol::resolve<cvector(const cvector&)>(Vector::cosh<complex>),
						sol::resolve<vector(const vector&)>(Vector::cosh<double>) )
		);
		module.set_function("pow", sol::overload(
						sol::resolve<double(double,int)>(std::pow),
						sol::resolve<double(double,double)>(std::pow),			
						[](const complex& base, int exponent) { return std::pow(base, static_cast<double>(exponent)); },
						sol::resolve<complex(const complex&,const double&)>(std::pow<double>),
						sol::resolve<cvector(const cvector&,double)>(Vector::pow<complex>),
						sol::resolve<vector(const vector&,double)>(Vector::pow<double>) )
		);
		module.set_function("exp", sol::overload(
						sol::resolve<double(double)>(std::exp),
						sol::resolve<double(double)>(std::exp),
						sol::resolve<complex(const complex&)>(std::exp<double>),
						sol::resolve<cvector(const cvector&)>(Vector::exp<complex>),
						sol::resolve<vector(const vector&)>(Vector::exp<double>) )
		);
		module.set_function("sgn", sol::overload(
						sol::resolve<vector(const vector&)>(Vector::vsgn<double>)
						,sol::resolve<double(double)>(Vector::sgn<double>) )
		);

		module.set_function("assoc_laguerre", [](int n,int m, double x){ return std::assoc_laguerre(n,m,x);});

		//module.set_function("interpolation", sol::resolve<double (const vector&,const vector&,double)>(Vector::interpolation<double>));

		module.set_function("interpolation", sol::overload(
											[](const vector& vec1, const vector& vec2 ,double x){
													return Vector::interpolation(vec1,vec2,x);
											},
											[](const vector& vec1, const cvector& vec2 ,double x){
													return Vector::interpolation(vec1,vec2,x);
											},
											[](const vector& vec1, const vector& vec2 ,const vector& x){
													return Vector::interpolation(vec1,vec2,x);
											},
											[](const vector& vec1, const cvector& vec2 ,const vector& x){
													return Vector::interpolation_complex(vec1,vec2,x);
											}
								)
		);

		module.set_function("fit_line", sol::overload(
											[](const vector& vec1, const vector& vec2){
													double a,b;
													Vector::fit_line(vec1,vec2,a,b);
													return std::make_tuple(a, b);
											}
								)
		);

		module.set_function("fit_plane", sol::overload(
											[](const vector& vec1, const vector& vec2, const vector& vec3){
													double a,b,c,d;
													Vector::fit_plane(vec1,vec2,vec3,a,b,c,d);
													return std::make_tuple(a, b,c,d);
											}
								)
		);


		module.set_function("fit_parabolic_line", sol::overload(
											[](const vector& vec1, const vector& vec2){
													double a,b,c;
													Vector::fit_parabolic_line(vec1,vec2,a,b,c);
													return std::make_tuple(a, b,c);
											}
								)
		);



		module.set_function("sort", sol::resolve<void (vector&)>(Vector::sort<double>));


		module.set_function("transform", sol::overload(
						[](const cvector& vec, sol::function fx){
							cvector res = vec;
							std::function<complex(complex)> stdfx = fx;
							for(auto& elm:res)
								elm = stdfx(elm);
							return res;
						  },
					   [](const vector& vec, sol::function fun){
								vector res = vec;
								for(auto& elm:res)
									elm = fun(elm);
								return res;
							  }
						 )
		);

		module.set_function("max", sol::overload(
						[](const vector& v){
							return Vector::max(v);
						}
				)
		);
		module.set_function("max_and_index", sol::overload(
						[](const vector& vec){
							double val;
							int idx;
							Vector::max_and_index(vec, val, idx);
							return std::make_tuple(val,idx+1);
						},[](const vector& vec,int nmin,int nmax){
							double val;
							int idx;
							Vector::max_and_index_with_range(vec, nmin-1, nmax, val, idx);
							return std::make_tuple(val,idx+1);
						}
				)
		);

		module.set_function("min", sol::overload(
						[](const vector& v){
							return Vector::min(v);
						}
				)
		);
		module.set_function("min_and_index", sol::overload(
						[](const vector& vec){
							double val;
							int idx;
							Vector::min_and_index(vec, val, idx);
							return std::make_tuple(val,idx+1);
						},[](const vector& vec,int nmin,int nmax){
							double val;
							int idx;
							Vector::min_and_index_with_range(vec, nmin-1, nmax, val, idx);
							return std::make_tuple(val,idx+1);
						}
				)
		);
		module.set_function("wrap_phase_deg", sol::overload(
						[](const double& v){
							auto c = std::polar<double>(1.0, v/180.*pi);
							return std::arg(c)*180.0/pi;
						}
				)
		);
		module.set_function("wrap_phase", sol::overload(
						[](const double& v){
							auto c = std::polar<double>(1.0, v);
							return std::arg(c);
						}
				)
		);

	}

} // namespace my_object

